Imaginary-time formulation of steady-state nonequilibrium in quantum dot models 
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We examine the recently proposed imaginary-time formulation for strongly correlated steady-state 
nonequilibrium for its range of validity and discuss significant improvements in the analytic contin- 
uation of the Matsubara voltage as well as the fermionic Matsubara frequency. The discretization 
error in the conventional Hirsch-Fye algorithm has been compensated in the Fourier transformation 
with reliable small frequency behavior of self-energy. Here we give detailed discussions for gener- 
alized spectral representation ansatz by including high order vertex corrections and its numerical 
analytic continuation procedures. The differential conductance calculations agree accurately with 
existing data from other nonequilibrium transport theories. It is verified that, at finite source-drain 
voltage, the Kondo resonance is destroyed at bias comparable to the Kondo temperature. Calculated 
coefficients in the scaling relation of the zero bias anomaly fall within the range of experimental 
estimates. 

PACS numbers: 73.63.Kv, 72.10.Bg, 72.10.Di 



Progress in nanoscale fabrication techniques has re- 
cently generated a great deal of interest in electron trans- 
port out of equilibrium. Well-controlled quantum dot 
(QD) structures in semiconductor devices have enabled 
thorough investigation of quantum many-body effects in 
confined geometry. One of the established phenomena is 
the zero bias anomaly (ZBA) due to the Kondo effect in 
Coulomb blockade devices^"— , where the scaling behav- 
ior of nonlinear conductance has been extensively veri- 
fied. Recently, more complex quantum dot systems such 
as molecular nano-j unctions^ and multi-channel Kondo 
dots^ with intricate device design have fueled intense re- 
search for electron transport mechanism in strongly cor- 
related regime. Nanoscale single-electron devices hold 
great promise not only in applications for quantum de- 
vices but also in development of general quantum many- 
body theory out of equilibrium. 

In the past few years, the strong correlation community 
has embraced the challenge of developing quantum many- 
body formulations out of equilibrium, which has lead to 
significant advances at various levels of theories. Unlike 
equilibrium quantum many-body theory where the ana- 
lytic and numerical theories play complementary roles, 
the nonequilibrium theory has only recently had full- 
fiedged numerical tools which could support or disprove 
the diagrammatic approximations, known as Keldysh 
technique^'^. So far, numerous algorithms have been pro- 
posed. However, most of the theories have yet to be fully 
established to have reliable predictive power in a wide 
range of strongly correlated transport. 

The main focus of the theoretical efforts has been 
the description of the transient behavior toward a 
nonequilibrium steady-state or the properties of an 
established steady-state. Here the steady-state con- 
cerns the dc current-carrying state driven by a time- 
independent bias in quantum dot geometry, as sketched 
in Fig. [TJ Quantum simulations of real-time be- 
havior of nonequilibrium steady-state have been per- 
formed using time-dependent numerical renormalization 
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FIG. 1: (a) Source and drain reservoirs for electron and a 
quantum dot level. Chemical potentials of the reservoirs are 
the same in equilibrium, (b) Voltage bias $ in steady-state 
nonequilibrium divided between chemical potentials of the 
reservoirs. The quantum dot energy level can be arbitrarily 
positioned with respect to the chemical potentials. 



group (INRG^'Si), time-dependent density-matrix RG 
(tDMRQi^iii), perturbative RG (PRQi^iH), functional 
renormalization group (fRG^^'^-), iterative summation 
of real-time path-integral method (ISPli^ii^), diagram- 
matic Monte Carlo ^^'^^ . Analytic methods of nonequi- 
librium Bethe ansata^S and perturbative steady-state ex- 
pansion^i have been developed. 

The imaginary-time formulation for steady-state 
nonequilibrium proposed by Han and Hearj*2£ takes quite 
a different approach from the above real-time techniques. 
Main motivation has been to extend the equilibrium 
quantum many-body theory within a similar mathe- 
matical framework of the imaginary-time formalism and 
therefore to easily adopt existing equilibrium numeri- 
cal techniques for complex quantum interactions such 
as molecular dots^'^. The method combines the equilib- 
rium many-body theory and steady-state nonequilibrium 
quantum statistics by extending the chemical potentials 
into complex variables, the Matsubara voltage. It has 
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been shown that, upon analytic continuation of complex 
chemical potential back to physical chemical potential 
the theory recovers the nonequilibrium dynamics, and 
the imaginary-time Green's functions map to real-time 
Green's functions. 

In this review, we show how the previously intro- 
duced spectral ansatz can be extended to include ver- 
tex corrections and give comprehensive discussions de- 
tailing the justification and the range of validity of 
the analytic continuation. We further enhance the 
computational method by modifying the conventional 
Hirsch-Fye^^ quantum Monte Carlo (QMC) algorithm by 
compensating the discretization error and obtain elec- 
tron self-energy with much improved reliability. These 
improvements lead to reliable differential conductance 
curves in accurate agreement with other existing meth- 
ods. We confirm that the Kondo resonance is destroyed 
at the bias of Kondo scale and the resulting scaling be- 
havior of conductance is consistent with experiments. 

The paper is organized as follows. In the follow- 
ing section HI we define the imaginary-time Hamilto- 
nian for nonequilibrium. We then show the equiva- 
lence of imaginary-time and real-time Green's functions 
in the perturbation expansion with arbitrary interaction 
through the analytic continuation. We give examples 
of different quantum dot geometry where this formal- 
ism can be extended. In the following section |lTl we 
briefly introduce the QMC procedure and propose an al- 
gorithm which fixes discretization errors in the nonequi- 
librium QMC. In section Hill we extend the previously 
introduced analytic continuation ansatz via Fade approx- 
imants and describe detailed numerical procedures to fit 
spectral functions. We make a direct comparison of con- 
ductance in the small to intermediate interaction regimes 
with other methods. In section HVl computational results 
are presented and compare them with existing theories 
and experiments. In the appendix, analytic structure of 
high order corrections to self-energy is discussed. 



I. IMAGINARY-TIME FORMALISM 

We start the discussion of nonequilibrium quantum 
theory with the understanding that the one of the funda- 
mental problems is that there are two different operators 
governing the quantum statistics and the time-evolution. 
In equilibrium, the time-evolution is given by e~**^ while 
the Boltzmann factor e"^^-^^^^' provides the quantum 
statistics with the chemical potential /i and the num- 
ber operator N. With a number-conserving Hamiltonian, 
[H, A^] = 0, the Boltzmann factor and the time-evolution 
operator commute. For this reason, the chemical po- 
tential /i is often set as the reference energy (/i — 0) 
without losing generality. However in nonequilibrium 
of quantum dot systems, bias voltage creates multiple 
electronic chemical potentials in the source and drain 
reservoirs. Once the tunneling between the quantum 
dot and the reservoirs is allowed and many-body inter- 



actions are turned on, the reservoir states mix with the 
time-evolution. Hence, the task of finding commutable 
nonequilibrium density matrix and the time-evolution 
operators becomes a challenging task. 

From the 1960s and 1970s, efforts have been made 
to formulate a Gibbsian statistical mechanics in steady- 
state nonequilibrium using the Liouville operator formal- 
ism'^^. In 1993, Hershfield^^ has revisited the problem 
in the context of the mesoscopic transport and has pro- 
vided a formal proof in a compact form of a density ma- 
trix for nonequilibrium steady-state. The idea consists 
of decomposing the full Hamiltonian in terms of the for- 
mal solution of scattering states of Lippman-Schwinger 
equation^!, and applying different chemical potentials to 
scattering states derived from each reservoirs. With the 
scattering state operator "^Ifccr "^ith the continuum in- 
dex k, reservoir index a = L,R (or , respectively) 
and the spin index ct, the nonequilibrium density matrix 
operator at the voltage bias $ is written as 



with the operator Y 

Y^^AkALka - V'k.V'fl.fc.). (2) 
ka 

However, despite its conceptual breakthroughs, Hersh- 
field's idea has not found practical implementations since 
fully interacting scattering states cannot be known a pri- 
ori. The method has been adapted only in the non- 
interacting models and in perturbative limits^*"'™. 

To overcome such difficulties, Han and Heary have 
proposed an imaginary-time formalism which constructs 
a formal perturbation expansion from a non-interacting 
nonequilibrium steady-state. The main advantage of the 
method over the original Hershfield's idea is that there 
is no need to construct the Hershfield's F-operator in 
the interacting limit. Many-body interactions are con- 
sidered rigorously in a perturbation series at arbitrary 
order built on the non-interacting nonequilibrium den- 
sity matrix g-0{Ho-'i>Yo) t^w^^ i^j^g yp-operator computed 
without many-body interactions. To compensate the dis- 
crepancy between the time-evolution and the Boltzmann 
factor, we introduce a mathematical trick of the Matsub- 
ara voltage. 

Electron tunneling of single quantum dot connected to 
source and drain reservoirs is modeled by 

^^0 = + Hqci + H^t (3) 

Hqc ~ ^ <^akC^ak<j^aki7 (4) 
ctk{7 

Had = ed^d^dcr (5) 

Hot = -^-^idlcaka + h.c.), (6) 

akcr 

with Hqc, Hod, Hot for decoupled reservoirs, quantum dot 
states, and QD-reservoir tunneling, respectively, c^^.^. are 
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the fermion operators for continuum state, rf^ quantum 
dot orbital operator, ta the tunnehng parameter and Vl 
the volume of the reservoirs. Then the non-interacting 
scattering state can be readily written down as^^i^ 



'k' 



n 



17] 



''a' k'<T 



(7) 



with the non-interacting retarded Green's function 
g°{u)) = [uj — €d + iT]^^ . r is the line- broadening of 
the quantum dot T = Tl Tj^ = 7r(t| + 4)iV(0) with 
the density of states of the reservoirs N{0). 

Even in the non-interacting limit, the problem of com- 
bining the density matrix and the time-evolution opera- 
tor remains. To resolve the issue, we extend the chemical 
potential into complex numbers such that the new chem- 
ical potential does not alter the quantum statistics while 
the time-evolution rate along the imaginary-time can be 
shifted. We introduce an effective non-interacting Hamil- 
tonian Kq, 



with the Matsubara voltage. 



for any integer m. 



(8) 



(9) 



With the addition of the Matsubara voltage, the density 
matrix given by Kq becomes e"'^^" = ^-i3[Ho+{iVra-t')Ya\ ^ 
Since [i/oj^o] — 0, the operator can be written as 
g-/3(ffo-*io) . ^-tPv^Yo^ ^j^j^ respect to the Fock basis 

constructed from the non- interacting operators ^f.^ , 
Hq and Yq are diagonal and the eigenvalues of Yq are half- 
or full-integers. Therefore, we have an operator identity 
^-iPv^Yo ^ I and 



-I3(Hq-^Yo) 



Po, 



(10) 



independent of i(pm- This effective Hamiltonian amounts 
to a electron system with the statistics given by po and 
the time-evolution is controlled at an independent phase 
velocity shifted by the pumping frequency of iipm- This 
extra complex term is only applied to the imaginary-time 
action and should not be interpreted as a physical damp- 
ing term in real-time formalism. 

Given the non-interacting Hamiltonian, Eq. ([8]) , we in- 
troduce the interacting Hamiltonian parametrized by 
iiPm ~ with the many-body interaction V as 

K = Ko + V = Ho + {t^m - $)fo + V, (11) 

and treat this as an effective Hamiltonian to the equilib- 
rium imaginary-time theory. In the following section, we 
will show that, upon the analytic continuation iip„i — > $ 
the thermal Green's functions map to nonequilibrium 
real-time Green's functions. 



The equivalence of the imaginary-time theory and the 
real-time Keldysh formalism can be best shown in explicit 
perturbative expansions. The real-time retarded Green's 
function for the quantum dot is defined as 

G"*(t) = -^e{t){{dH{t),d^H{0)}) (12) 
= e{t)[G>{t)-G<{t)], (13) 

with the lesser and greater Green's functions defined as 



G>{t) = -i{dH{t)d^Hm 

G<{t) = z(4(0)d^f(i)). 



(14) 
(15) 



The subscript H refers to the evolution in the Heisenberg 
picture with the full Hamiltonian, d''{t) — e**^(i^e~**^. 
6{t) is the step-function. 

Usually an interaction picture is defined with a time- 
dependent many-body interaction V turned on adiabat- 
ically from the infinite past {t — T as T — > — oo, see 
Fig. [5]). The Green's functions can be defined in the 
Keldysh contour as 



G>{t) 



-iZ, 



0,neq 



Tr 



(16) 

The real-time evolution is given in the interaction picture 



dHt) 



-itHo 



and Vi{t)e'*"'>Ve 



-itHo 



Time 



variables are defined on the Keldysh contour, denoted as 
K. The nonequilibrium partition function is defined as 

Zo,neq — Trpo- 

Here, to avoid explicit time-dependence to the interac- 
tion V, we use Gell-Mann and Goldbcrgcr's'^*' construc- 
tion of steady-state when the limit T — > —oo is taken. 
Their formalism uses full-strength interaction instead of 
adiabatic interaction with the turn-on time T integrated 
from the remote past to the present. For example, a 
scattering state operator defined as 



4^ 



fe' 



satisfies the Lippman-Schwinger equation 



4 



1 



ek- C + ir] 



[V + Hot,cl], 



(17) 



(18) 



with the tunneling part of the non-interacting Hamilto- 
nian Hot- Here the infinitesimal 77 determines the di- 
rection of the time propagation. The Liouville opera- 
tors are defined as £A = [H,A], CqcA = [Hqc,A], and 
g-j£T^ = e~^^^ Ae^^'^ , etc. Here the adiabatic limit 
T — > —00 is replaced by an integral via rj J_^dT e^'^ 
where the start time T of interaction is averaged over 
the whole range (— oo,0]. This averaging effectively can- 
cels out the transient phase of each scattered wave which 
has propagated with the interaction turned on at time T. 

The greater Green's function G^{t) can be represented 
as Fig. [5] with a d-electron created at time on the up- 
per contour and an electron destroyed at time t on the 
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{a) 
t^T 



(6) 



FIG. 2: (a) Keldysh contour for greater Green's function with 
the first order scattering (marked by X) happening between 
time t = and t [corresponding to the case (i) considered in 
the text] (b) Time ordering of case (ii). (c) Case (iii). 



lower contour. Then, the first order perturbation correc- 
tion from Eq. (|16p can be decomposed into three different 
time-ordered terms according to the time of interaction 
t' as (i) 0+ < t' < t_ (ii) T+ < t' < 0+ (iii) t- <t' <T_, 
as shown in Figs. [lUa)-(c), respectively. The contribu- 
tion for (i) can be expressed as (with Z^^^^ omitted for 
brevity) 

"'0 „„,,, 



^it{E„-E^) 



t{E„-Ek) 



mfc"fc„ 



nmk 



E„ 



m 



where the states denoted by n, m, k are Fock states 
constructed from the non-interacting scattering states. 
(• • • )o is defined as Tr[po • • • ]• Note that zero of the en- 
ergy denominator Em — Ek does not lead to a singularity 
since e'*^^"-^'") - e»*(^"-^'=) = 0. 

For the perturbation occurring in the interval extend- 
ing to the infinity on the upper time contour [Fig. Elb)] , 
the integral for (ii) becomes 



-IT] 



dTe^'^ [ dt' {e- 
Jt 



tHr 



de 



-itH, 



itAE„,„ d.nmdl^j,Vkn 



(ii) and (iii) combine to 



itAE„ 



Po,me 



„itAE„ 



nmk 



AE„ 



(22) 



Finally the first-order perturbation to the greater Green's 
function becomes 



E 



PO.,' 



JtAE„ 



JtAE„ 



dnrrS^rnkd 



AE, 



mk 



E [PCne^'^^"-^- -Po.^e^*^^-^] ^"-^-fe^fc» (.23) 



AEr, 



irj 



We perform the same perturbation theory to the ther- 
mal Green's function defined with the imaginary-time un- 
der the Hamiltonian K, Eq. (|Tl1) . as 



(24) 



g{t) = -(r.dK(r)4(o)), 

with the time-evolution given as dx (t) — 
ing the interaction picture, the Green's function is ex- 
panded in a perturbation series as 



Us- 



Git) 



Po 



r,e-^o'^Hr')dr'^(^)^t(o)l ^ (25) 



where the time-evolution in the interaction picture is 
given as d(r) = e'^^°de~'^^° . For r > 0, the first per- 
turbation contribution due to the scattering of V at r' 
can be grouped into two cases; (i) < r' < r and (ii) 
T < r' < /3. 

The contribution from (i) < t' < t is 



Tr 



= E ^0'" t^'' 



e-(^<-r)/forfg-(r-r')Koye-r'K„^t 
dnm^mkdi^^ 



AK„ 



..tAK„ 



AEh„ 



irj 



/o 
(20) 



nmk 



AK„ 



dr' 
(26) 



Po- 



with AEnm = En — E„i ■ Similarly for the case (iii) with 
the interaction on the lower contour, we have 



-IT] 



with the eigenvalues Kn of Kq and AKn 
Here we have used the key relation Eq. (flUl) . e' 
If the analytic continuation iipm — > $ is formally carried 
out followed by T -> it, Kn En and the above integral 
for < t' < r becomes identical to the real-time Green's 
function Eq. ((TO)) for 0+ < < Now, for the case of 
(ii) T < t' < (3, the integral becomes 



/3 



Tr 



E 

nmk 



AEn 



irj 



(21) 



EL 

nmk 



Po,ne 



tAK„ 



Po,me 



tAK„ 



AKnm 



dr' 
,(27) 



By denoting the time ordering on the Keldysh con- 
tour as {ab)K for an event a following b, the Eqs. ([T9]- 
can be represented by a cyclic permutation of 
{{dVd^)K,{dd^V)K,{Vdd^)x}, respectively. Rearrang- 
ing the indices in Eqs. (|20II21I) . the contributions from 



which transforms to the real-time result, Eq. ([2T|) . with 
the analytic continuation except for the adiabatic factor 
irj. 

We discuss the subtlety of the adiabatic factor irj, the 
presence of which is only relevant for the energy shell of 



(a) \n) 




(a) 



(b) 



id) 



FIG. 3: (a) Interaction V mapping a state |n) to |m). (b) 
Initial and final states without changes in the Vo-number of 
scattering states, (c) Initial and final states with different 
Yo-numbers. (d) Contraction with the d-creation and anni- 
hilation operators. Each line represents the contraction of 
scattering state basis (V'lfeCTV'afco-)o and contributes the fac- 
tor t^\g{€ak)\'^ ■ Two diagrams are with |n) and |m) states 
interchanged. 



En = Em in Eq. 



PO,: 



itAE„ 



nmk 



(28) 

In equilibrium, po is only given by energy and the above 
expression is zero since po,n — Po,m = 0. To extend the 
argument to nonequilibrium, we take an explicit example 
of one-quantum dot with a two-body interaction such as 
the on-site Coulomb interaction V = Uud^ndi- In terms 
of the exphcit scattering state operators Eq. ([7]), 



4 

V 



{a,k} 



iaiff*(ei) 



ia35*(e3)^t 



a3k3-l 



(29) 

V'a4fe4i ) • (30) 



As depicted in the Fig. El^b) , if the incoming and out- 
going states |n) and \m) conserve the eigenvalues for lo 
operator (Yo-number), the factor (po,,i — po,m) is zero. 
Figure ^c) shows the terms in V [Eq. (|30l) ] which do 
not conserve the lo-number. However, when expecta- 
tion values are taken between the states \n) and |m), it 
can be shown that there is a counter-term which leads 
to a cancellation. As demonstrated in Fig.[3Jd), the con- 
nected legs represent the contraction {ijjl^^^ipaka)o and 
contribute the factor t'^\g{eak)\'^ from Eqs. ^ and (|30l) . 
If the source and drain reservoirs are given by the same 
continuum density of states, there is always a contribu- 
tion with the same magnitude from when the states \n) 
and |m) are interchanged, therefore leading to the can- 
cellation of the factor {po,n — Po,m)- We emphasize that 
the above argument does not require assumptions for the 



(c) 



(d) 



FIG. 4: Nonlinear Transport through (a) single QD, (b) side- 
coupled QDs and (c) parallel QDs can be studied within the 
imaginary-time formalism, (d) Green's functions in serially 
coupled QDs may not be correctly continued to lesser Green's 
functions. 



source-drain symmetry — tji or the particle-hole sym- 
metry. It also holds for different types of on-site interac- 
tion. 

To summarize, the time-ordering of the real- 
time greater Green's function can be matched to 
the imaginary-time-ordering [denoted by (• • • )/] as 
{dVd'')K ^ {dVS)i and {Vdd'^)K + {ddW)K ^ 
{ydd'^)i. The higher order contributions can be checked 
in the similar manner to the first order. For exam- 
ple, in the second order, the time-orderings in the real- 
time theory can be matched as {dVVd'^)K ^ {dVVd''')j, 
{VdVd^)K + idVd^V)K ^ iVdVd^)i, and iVVdd'')K + 
{VddWjK + {ddWV)K ^ {VVdd'f)!. Such topologi- 
cally equivalent graphs between the imaginary-time and 
real-time expansions at each perturbation order are ex- 
pected since the two theories are known to be equivalent 
in equilibrium. 

The lesser Green's function can be shown to be equiv- 
alent to the imaginary-time Green's function G{t) for 
T < 0. As shown previously^, if Fourier transforma- 
tion is performed on the real- and imaginary-time Green's 
functions, the retarded Green's function G"''^'(w) can be 
obtained by analytically continuing the thermal Green's 
function fj(za;„) via za;„ ^ td + i?]. 

So far, we discussed the analytic continuation in single- 
QD structures [Fig. [^a)]. However, the formulation can 
be straightforwardly extended to much wider range of 
quantum dot structures such as the side-coupled QD 
and parallel-coupled QD systems as shown in Fig.|3l^b-c). 
For serially coupled QD systems [Fig. UJ^d)], the current 
imaginary-time formulation does not have an analytical 
continuation to real-time Green's functions and hence the 
formulation in this work cannot be applied. 



II. COMPUTATIONAL ALGORITHM 

We implement the imaginary-time formalism for a nu- 
merical nonequilibrium technique using quantum Monte 
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Carlo method. We use the Hirsch-Fye (HF) algorithm-- 
which uses only the interacting orbitals (QD site) as 
the dynamic variable after integrating out the non- 
interacting reservoir electronic states. The QD Green's 
function is stochastically updated along the discretized 
imaginary-time lattice using the local update method 
by Blankenbecler et al^. Calculation of nonequilibrium 
interacting Green's function does not require any main 
modifications of the standard (equilibrium) QMC code. 



In this work, we analytically continue the retarded 
(self-energy) functions instead of the lesser/greater 
Green's functions as functions of frequency since the re- 
tarded functions have simpler analytic structure of ratio- 
nal functions, as opposed to the exponential functions for 
the lesser/greater Green's functions as functions of time. 

B. Quantum Monte-Carlo method and Self-energy 



A. Non-interacting Green's function 

Independent sets of simulations are performed with dif- 
ferent iifm^s in Eq. treated as fixed parameters to 
QMC. In the following calculations, < m < 5 have 
been used. Given iipm, the non- interacting QD Green's 
function can be easily derived as 



1 



ak 



iuJn — Ko^ak 
l5(ea/c)P 



ak 



lUJn - eak- a 



iV ■ Sign(Im2;„m) ' 



{'4^cxka\da) (31) 

(32) 
(33) 



2 



with Znm 

of iipm ^ followed by iw„ — > lo + irj is performed, 
the thermal Green's function in the Matsubara frequency 
transforms to the real-time retarded Green's function. 



E 



1 



uj + ii] + ?T uj + irj + iV 



(34) 



Fourier transformation to the imaginary-time variable 
(r > 0) gives 



(35) 



E ||g°(eafe)pe--[^-+"('^"-*)/2l[l - Um 



with the Fermi-Dirac function in the a-reservoir, /«(£) = 
[1 -I- e^'-'^""*/^-']"^. This expression is later used as the 
input to the QMC calculation. If we perform the analytic 
continuation i(pm — > $ followed by r — li, the Green's 
function in the imaginary-time transforms to 

- E Sl5"(^afc)Pe-'*^-[l - Ue^,)], (37) 



ak 



which is nothing but the non-interacting greater Green's 
function of QD orbital, G^{t), in the steady-state 
nonequilibrium. In this sense, the thermal Green's func- 
tion before analytic continuations contains full informa- 
tion of the retarded, greater and lesser Green's functions 
in the real-time formulation. 



The QMC procedure follows the standard Hirsch-Fye 
algorithm with the initial Green's function Eq. ([55]) at 
fixed i(pm- QMC method stochastically samples the 
fcrmionic phase space via auxiliary fields in the action in- 
troduced by the Hubbard-Stratonovich transformatior^. 
The auxiliary fields are updated according to the effective 
Boltzmann factor—. The only modification to the HF al- 
gorithm is that the Monte Carlo (MC) Green's function 
G(r, r') at an auxiliary field configuration is complex in 
contrast to the equilibrium calculations. Since G{t,t') 
is complex, the ratio of Boltzmann factors for the new 
and old auxiliary field configurations is also complex in 
general. Therefore for any observable (A) we compute 
the ensemble average as 



(A) 



E„/(n)A(n) _ j:n^"'"Mn)\f{n)\ {{e^' A)) 



aie\ 



(38) 

with the effective Boltzmann factor f{n) for a auxil- 
iary field configuration n and its phase factor e*^" = 
f{n)/\f{n)\. The ensemble average ((• • • )) is taken over 
the Markov chain of the Monte Carlo configurations cho- 
sen by the probability |/(n)|. In the calculations shown 
later, the average phase factor |((e*^))| remained close 
to one (typically 0.9 — 1.0) and the statistics has been 
quite robust. We note that at $ = and iipm — 0, the 
QMC calculation is completely identical to the equilib- 
rium QMC method. 

The QMC Green's function defined on a discrete 
imaginary-time mesh = iAT(AT = (3/N, i = 
0, - ■ ■ ,N — 1) is updated by the QMC Dyson's equa- 
tioi*2i^ 



G = Go + (Go-/)(e^-/)G 



(39) 



for the Green's function matrix defined as G^ = 
G{Ti,Tj). Here V represents the auxiliary-field coupling 
to electrons after the Hubbard-Stratonovich transforma- 
tion^'*''^^ of many-body interaction. The diagonal com- 
ponent of the matrix G is chosen as the greater Green's 
function by'^-^ 



(40) 



This choice of zero-time Green's function introduces dis- 
cretization error when the time variables are integrated 
in the above Dyson's equation or in other observables. 
In nonequilibrium calculations, the Green's function, 
Eq. (|36p. has an additional oscillation due to the iifm 
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dependence, and the systematic error of discretization 
becomes more significant. Spurious structures in the low 
frequency self-energy observed in Re^ are attributed to 
the discretization error^^, which can be confirmed in com- 
parison with the continuous-time QMG^^''^^. In Figs.[5][a- 
b), the discretization error is compared for At = 1/5 and 
1/10. The energy unit is the non-interacting broadening 
r. Although the self-energy at high frequency uj„ is well 
convergent, the low frequency E(ia;„, shows discon- 
tinuous jumps at a;„ « as ipm increases. This kink be- 
comes smoother as Ar becomes small. Without the cor- 
rection, the analytic continuation misinterprets the kink 
in the self-energy due to incoherent spectra and exagger- 
ated the destruction of Kondo resonance at finite bias^^. 

By adopting a similar trick for Fourier transforma- 
tion considered in the continuous-time QMC''^'''®, the 
Green's functions in the discrete-time QMC has been 
measured as follows. When involved in a time-integral, 
we use the non-interacting Green's function matrix Gq 
with the diagonal elements augmented by Go a = 
\ [Go(t, +0+,T,) + Go(t, - 0+,r,)] = Go(T, + 0+,T;)-i, 
i.e. Go = Go — \l- Then the Dyson's equation can be 
rewritten as 

G = Go + GoSGo + i (sGo - GqS) , (41) 

with the S-matrix defined as 

S = (e^ - /)GG(7^ (42) 

Through Monte Garlo updates we measure the Fourier 
transformed S-matrix as 




(43) 

Here the matrix G is calculated at each update of the 
auxiliary fields and G^^^ is calculated and stored at the 
beginning of computation. The last term in Eq. (|4ip 
vanishes due to the time translational symmetry and 

Q{ibJn) = Q^(ibj,^) Qa(ibJn)S(ibJn)Qa{i'^n)- (44) 

As shown in Fig. Eljc) , the unphysical structure at w„ « 
disappeared even at At = 1/5 after the discretization er- 
rors have been corrected. Figure [5jd) at a finite bias $ 
shows less curvature in the self-energy, which suggests 
that the correlation effects become weaker as bias in- 
creases. Green's functions evaluated this way showed ex- 
cellent agreement with the continuous-time QMG^ at 
low fermion frequencies ia;„ at large Matsubara frequen- 
cies iifra with computationally accessible At. In the cal- 
culations presented below used At = 1/5 unless men- 
tioned otherwise. 

The self-energy of the QD Green's function 
'Ej{iujn,i^m) is then computed in the same manner 
as in equilibrium theory, via the Dyson's equation 

T.{iLOn,ifm.) = [Q^{i>^n)\ ^ -~ [^^(«W„)]"^ . (45) 



1 ° 



(a) No correction W\\^;^3^s^^^g 
* = 0, At=1/5 1 


(b) No correction M\\^^^~S^^^^^0S^ 

4 = 0, At= 1/10 


(c) With correction Vv;i|^^^p>>^ 
* = 0, At=1/5 VNi*"**^ 


(d) Witt! correction ^^^^g^^^^ 
- «=1,At=1/5 



-8 -6 -4 -2 2 4 6 8-8 -6 -4 -2 2 4 6 S 
iVIotsuboro freq, tvlotsuboro 'req, uj^ 



FIG. 5: Imaginary-time electron self-energy at [/ = 10 
and /3 = 24 for Matsubara voltages ifrn with m = 
(filled circle), • ■ • , 5 (open circle), (a) Conventional Hirsch- 
Fye algorithm without the discretization correction shows 
spurious structures at small uj„ and finite ipm- (b) With 
smaller discretization At = 1/10, the spurious structures 
become weaker, (c) The discretization correction produced 
smooth self-energy at small uj„. (d) At higher bias $ = 1, the 
curvature at high ipm becomes weaker, suggesting reduced 
correlation efi'ects. The unit of energy is given by the non- 
interacting level width of QD, F = 1. 



This self-energies I](iw„,«<y9m) at different iipm values are 
computed in separate sets of QMG runs at each iipm- 
The numerical data for T,(iLUn,i(pm) is analytically con- 
tinued to the retarded self-energy T,^'^*{uj) with the real- 
frequency ui. The numerical procedure will be fully dis- 
cussed in the next section. 

We comment on why we choose to analytically continue 
the self-energy instead of the Green's function directly. 
As will be clear in subsequent discussions, the analytic 
form of the perturbative energy self-energy is more read- 
ily written down, which makes the numerical procedure 
more transparent. From the numerical standpoint of the 
discrete-time QMG, the discretization makes the high 
Matsubara frequency data less reliable for T,{iuJn,i(pm) 
and Q(iujn,iLPm)- However, the systematic errors in an 
analytically continued S(w) at large lo are less problem- 
atic since the frequency term uj in the Dyson's equation 
G''^*(a;) = - Erf -f iF - T.{uj)]-^ dominates S(a;). 



III. NUMERICAL ANALYTIC CONTINUATION 
A. The spectral ansatz 

It seems a formidable task to perform an analytic con- 
tinuation on numerical data in I](iw„, i(p,„). To guide 
the analytic continuation to a correct form we start with 
the self-energy in the second order of Goulomb interac- 
tion in the Anderson model, with the diagram depicted 
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ei , u„ — v„ 



FIG. 6: Self-energy diagram of second order perturbation in 
the Coulomb parameter I] of the Anderson model. 



in Fig. |S1 The retarded self-energy can be easily calcu- 
lated from Y^(t) = U\G^(t)fG^{-t) and S'-^*(t) = 
Q{t)'^>{t) - E<(i)] with the step-function Q{t), 



^ r [ 



^ /l(l-/2)/3 + (l-/l)/2(l-/3) ^ 

w -I- iry - £1 + £2 - £3 

with the shorthand notation /j = fai{£i)- Po(e) is the 
non-interacting QD spectral function. If we do the same 
diagram in the imaginary-time formalism with the Hamil- 
tonian K, the self-energy is 



r f 



^, /l(l-/2)/3 + (l-/l)/2(l-/3j^ ^^ 

iujn - £1 + £2 - £3 

with ei = Ci -I- ai{i(pm — *i')/2. Here we have used the 
relation for the Fermi-Dirac function, 

f{e + a'-^)=f{e-a^)^Ue), (48) 

which is equivalent to Eq. (jlOp . By combining the reser- 
voirs indices 



7 = ai — q;2 + as 



(49) 



and the energy of an electron dressed by an electron- 
hole pair as £ = £i — £2 + £3, we can rewrite the above 
expression as 



E 



7=±1,±3 



de- 



with the spectral function cr^(e) defined as 



(50) 



Ql — Q2+Q:3— 7 



^ r f 



x[/l(l-/2)/3 + (l-/l)/2(l-/3)] 

x<5(a; - £1 + £2 - £3). (51) 



In the second order of interaction the self-energy spectral 
function aj{e) is independent of iipm ~ ^- However in the 
higher order of perturbation, it is no longer the case and 
we need to incorporate the ii^m — ^ dependence in the 
spectral function as 



E(iw„,i(p,„) = 



d£ 



a^{e)Q^{t, iifim - <&) 



$) 



-l=oddZ'' """" 2" 

(52) 

with any odd integer 7. Q-y(£, i^Pm — ^) is the correction 
due to higher order diagrams. See the appendix for de- 
tailed discussions for this generalization of the spectral 
representation and its analytic properties. We approxi- 
mate the function Q-y by a Fade approximant 



y(£,z) 



(2), 



l + D''^\e)z + D\'>{e)z-^ + 



(2), 



(53) 



Therefore we seek the best spectral representa- 
tion of the QMC-computed self-energy by treating 

{(j~j{e),C^\e),D'i^\e)} as fitting parameters. In the 
following calculations, we limit the 7-branches to 7 = 
±1,±3,±5,±7 and the Fade approximants to the first 
order n = 1, which already required fitting 24 functions 
simultaneously. We will discuss in the next section the 
effects of the Fade approximants. We emphasize that, al- 
though the Fade coefficient functions {C^\e) , D'"-y^\e)} 
are adjusted in the numerical fit, they do not contribute 
to the (real-frequency) self-energy after the analytic con- 
tinuation iLp,n — > 



ImE'-'=*(w) = -7r^CT^(w). 



(54) 



Real part of (w) is obtained from the Kramers-Kronig 
relation. 

The electron self-energy satisfies the general relation 



E(iw„,i((£'„ 



P(- 



-i</5m)]* 



(55) 



as can be seen in the non-interacting Green's function, 
Eq. p2p . For a particle- hole symmetric system, the non- 
interacting Green's function Qgiiuinjifm) in Eq. p2p is 
invariant with iipm — —iipm + $ and we derive sym- 
metry relations 



(56) 



In the previous work , the relation a^{e) — a f{~e) has 

been incorrectly applied as a^{e) = a ^{e) and this led 

to overly rapid reduction of the Kondo resonance at finite 
bias. As pointed out^"^ later, correct constraint produced 
a good agreement with other method^ in the moderately 
interacting limit, U = 5 [see Fig. 2(d) in Ref.— ]. In this 
work, we do not impose any symmetry relations in the fit, 
and the resulting spectral functions recovered the above 
relations numerically. 
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B. Fitting procedures 



Calculation of conductance 



With the above spectral ansatz, we perform the least- 
square fit to the numerical self-energy generated by the 
QMC calculations with defined as 



JV-l 



M 



nm 



n=-N m=0 



'>^QMc{i^miVm) 



(57) 



The deviation between the self-energy fit T,fit{iuJn,i(pm) 
in the above ansatz and the QMC generated data 

TjQMc{i'^n,i'fm)- We fit more accurately the low fre- 
quency self-energy data with an effective cutoff function 



.-2 



The normalization factor A/" is defined as 



N-l M 

E E 

n—-N m— 



(58) 



(59) 



In the following calculations we used M = 5 and lon = 

8.or. 

In QMC applications, the analytic continuation has 
been one of the most controversial topics. Since the 
transformation of imaginary-time data to real-time infor- 
mation is an ill-defined procedure, noise in the imaginary- 
time data can lead to severe uncertainties in spectral 
functions. In the past, several algorithms have been pro- 
posed and the maximum entropy method based on the 
Bayesian inference^ and the method of stochastic image 
generation^ have been widely used. In this work, we 
have not incorporated such methods where the focus so 
far has been limited to finding right spectral represen- 
tations. Simultaneously finding a fit to many functions 
(24 as noted above) has already been quite an extensive 
task computationally. Refining the analytic continuation 
method remains an important area of future improve- 
ment for the imaginary-time theory of nonequilibrium. 

We discretize the frequency on logarithmic mesh sys- 
tems with 201 frequency points centered aX u — Q 
over the energy range [—30,30]. Minimization of is 
achieved iteratively by using the Newton's steepest gra- 
dient method^^ . Once the fit reaches a certain threshold 
of accuracy {\f^ < 0.06), we regularized the spectral 
functions through third-order polynomial smoothing to 
reduced unwanted noise. The smoothing has had mostly 
insignificant effects and often has been unnecessary. The 
fractional error \/)(^ in the fit resulted in the range of 
1 — 6 %. Generally, the Pade approximant term pro- 
duced better fits. Due to the dense frequency points near 
uj — 0, the update of spectral functions at small frequen- 
cies tends to be very slow. Therefore, for faster conver- 
gence, we used adjustable mesh systems which evolved 
from a coarse to a dense frequency mesh as the iterations 
progressed. 



Once the retarded Green's function is obtained, the 
current in a single-quantum model can be computed from 
the Meir-Wingreen's formula^. 



1=1^1 &A(e,$)[/L(e)-/«(6)], 

with the QD spectral function at the bias $, 

1 



(60) 



(61) 



The differential conductance G is obtained from numer- 
ical differentiation of the current by 



dl 



(62) 



The differentiation has to be taken on discrete values of 
$i (i = 0, 1, • • • ). For $i = 0, the linear conductance is 
obtained from 



G($==0)/Go-7rr / deA{e,0) 



dm 

de 

with the conductance quantum Go given by 

2e2 



Go — 



h ' 



(63) 



(64) 



For the first non-zero bias (i = 1), we evaluate the deriva- 
tives from a third-order polynomial at $ = $i deter- 
mined from the values of current {— /2, — /i, 0, /i, at 
bias {-$2,-$i,0, $i,$2} as. 



3 A$ 

with $i = iA$ (i = 0, 1, 2). For higher bias <!>,; {i > 1) 



(65) 



G. = i 
2 



li+l — li 



h-1 - h 



(66) 



D. Comparison to other methods 



To demonstrate the validity of the imaginary-time 
QMC (ITQMC), the differential conductance (G = 
dl/dV) is compared to other methods where data is 
available. Comparison to other methods of fR( ] | i4^i5,i7 ^ 
ISPli^ilI, tDMRG" is shown for the weakly interact- 
ing limit in Fig. [7{a). In such limit the spectral ansatz, 
Eq. (|50p . becomes an exact representation of the self- 
energy at all bias and the resulting conductance is in 
good agreement with other methods. Even in the inter- 
mediate coupling limit U/T = 5 in (b), the comparison 
to the time-dependent numerical renormalization group 
(tNRG^) results is very good. Based on this concrete 
comparison and numerical efficiency of QMC for tack- 
ling complex QD models in strongly correlated regime, 
this imaginary-time method provides an efficient tool in 
nonequilibrium transport theory. 
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ITQMC 






fRG 






ISPI 




• 


tDMRG 








(a) U/r= 5 







-•-ITQMC 

o tNRG 



1 .4 



(b) u/r= b 



Bias, */r 



Bias, $/r 



FIG. 7: (a) Comparison of difFerential conductance in 
weakly interacting limit ai U /T — 3 from the imaginary-time 
QMC (ITQMC) at P — 25, functional renormalization group 
(fR G^^'^^'^^ , data taken from''^^), iterative summation of path 
integral (ISP I^^'^^ ), time-dependent density matrix renormal- 
ization group (tDMRG(^). (b) Conductance in the interme- 
diate coupling regime at U/T = 5. Curves are from ITQMC 
{/3 — 24) and the time-dependent numerical renormalization 
group (tNRG^ at 13 = 25). 



IV. RESULTS AND DISCUSSIONS 

The spectral functions of the self-energy from the 
ansatz Eqs. (|52II53I) are shown in Fig. [5] and they demon- 
strate the nature of quasi-particles dressed with particle- 
hole pairs. The parameters are U = 10, /3 = 36 and 
<f> = 0.5. In addition to the Kondo resonance at zero 
frequency, there are structures in the spectral functions 
at w = $/2 for 7 = 1, cj = 3$/2 for 7 = 3, etc. For 
7 = 1, the sum of reservoir indices of a dressed electron 
has J2i ^ — \ ^"^^ this leads to a resonant structure 
at the electron frequency measured from the combined 
chemical potential i$. Similarly, there is a resonant 
structure at w = |$ for 7 = 3. These effects of cross-lead 
particle-hole excitations appear as shoulders in the spec- 
tral function A{u)) in Fig. 13 The magnitude of the Fade 



terms c'^^ (w) and D'^''' (uj) is much smaller than one, and 
this suggests that the higher-order Fade approximants 
will not significantly change the data presented here. We 
also note that the symmetry relations Eq. (j56p have been 
numerically verified. 

Spectral functions in the strongly correlated limit U — 
10 and 14 are shown in Fig. [SI In (a), the Kondo reso- 
nance develops sharply on top of the incoherent charge 
excitations with the Hubbard peaks at a; ~ ±C//2. The 
inset shows the spectral function for the whole frequency 
range. The bias values are $ = 0, 0.2, 0.6, 1.0, 1.5 and 2.0 
(top to bottom curves). For clarity the curves are shifted 
vertically by 0.05. The Kondo resonance is strongly 
quenched as the bias is applied. There appear spec- 
tral shoulders at w = $/2 (single vertical lines) and at 
ijj = 3(f>/2 (double vertical lines). Their strength is con- 
siderably weaker than reported in other work a^'^^ . 

The Kondo temperature is estimated from the half- 
width-at-half-maximum (HWHM) of the Kondo peak at 
zero bias $ = 0. Due to the incoherent charge back- 
ground at T:TA{bj) w 0.2, we read off the HWHM at 



1(1) 




(b) Ee. = V2 



(c) Ee. = 3V2 



FIG. 8: (a) Spectral functions for the ansatz Eqs. (|52I53I) at 
U = 10, /? = 36 and $ = 0.5. Dashed lines are for negative 

7's, eg. (T_i(tj), o"_3(ij;), etc. As expected, o"^(e) = o" y( — e) 

is satisfied in the particle- hole symmetric limit. Short vertical 
lines indicate cj = "I'/2 (single) and ui = 3<l>/2 (double line), 
(b) Intra-lead single-particle excitations dressed by particle- 
hole excitations. They contribute to (Ji{uj) at excitation en- 
ergy of "l>/2. (c) Inter-lead excitations for 173(0;) at excitation 
energy of 3'l>/2. 



tiTA{ujk) =0.6 and ojk = 0.075. The Kondo tempera- 
ture from the renormalization group (RG) theory in the 
strong coupling regime^i^ has 



-K 



with the HWHM 



exp 



ttC/ 
ST 



7rr\ 
2U) ' 



at 



,RG 



= -T, 



RG 



K 



= 0.066, 



(67) 



(68) 



in a reasonable agreement with our numerical estimate. 

Nonlinear conductance for [/ = 10 is shown in Fig. [10] 
(a) with Fade approximants to the first order and (b) 
without the Fade correction. The comparison demon- 
strates that the corrections are insignificant, at least in 
the particle-hole symmetric Anderson model. Since the 
current is an integral of the spectral function, details in 
the spectral weight shift tend to be insensitive to differ- 
ent approximations of analytic continuation. The broken 



11 




-2-10 1 
frequency, w 



FIG. 9: (a) Spectral functions of one-particle Green's function 
at [/ = 10, /? = 36 and bias voltage $ = 0, 0.2, 0.6, 1.0, 1.6, 2.0. 
The Kondo resonance becomes quenched quickly as the bias 
$ is applied. The spectral functions are shifted by 0.05 for 
clarity. Short vertical lines mark spectral features at oj = $/2. 
Inset: spectral function for the whole frequency range, (b) 
Spectral functions at i7 = 14, /3 = 48. Double vertical lines 
denote the spectral contributions at uu = 34'/2. At larger U 
the Kondo peaks get quenched faster to lower intensity. 



lines in (a) are derived from Eq. (|60|) with the equilibrium 
spectral fmiction Aeq{Li>) calculated at $ = 0. Therefore, 
the reduced ZBA width (about 50 - 60 %) in the full 
nonequilibrium calculations are due to the destruction of 
the Kondo resonance at finite bias. 

For a comparison to experiments and other theories, we 
estimate first the temperature T1/2 at zero bias at which 



the linear conductance becomes the half conductance 
quantum, G(Ti/2,$ = 0) = ^Gq. T^ji « 1/14 = 0.071 
at a similar energy scale with the above HWHM ujk 



0.075 and 



,RG _ 

■'k — 



0.066. Then, from Fig. ilja), we esti- 
mate the bias for half conductance quantum at the min- 
imum temperature of simulation T,nin, G(T,nin, $1/2) = 
\Gq is estimated to be $1/2 ~ 0.135 at Tmin = 1/60. 
From the phenomenological scaling forn*^"^ of the con- 
ductance in the leading order of temperature and bias is 
written as 



G(T, $) 
Go 



= 1 



Ct 



T 



- cv 



(69) 



By definition, ct = 1/2 and cv can be derived by solving 
G(Tinin = 1/60, $1/2) = 5G0. Then we have an estimate 
for the ratio of the coefficients a at C/ = 10 as 



-'1/2 



($ 



1/2 



0.26. 



(70) 



Similar calculations have been repeated for U — Yl and 
14 and the results are summarized in TABLE HI Due to 
the small Kondo temperatures at large ?7, the maximum 
linear conductance at zero bias reached short of the con- 
ductance quantum at G/Gq = 0.73 for U = 12,/3 = 60 
and G/Go = 0.63 for = 14,/? = 60. The QMC esti- 
mates for Ofgnic are about (Xqmc ~ 0.2. We note that the 
incoherent spectral background in the Anderson model is 
at 0.1 — 0.2 in Fig.[9]with the conductance also approach- 
ing 0.1 — O.2G0 at high bias in Fig.[lUl as opposed to theo- 
retical predictions from the Kondo model or renormalized 
resonant level model. Therefore the above estimates of 
T112 and $1/2, hence (Xqmc, should be taken with some 
caution when compared to other theoretical models. We 
also note that the estimates of cnqmc have been derived 
from finite values of T1/2 and $1/2, instead of taking the 
limit T, $ -)> OM. 

For large /7, the QMC calculations tend to produce 
overestimated Kondo resonance HWHM, 10 Ki compared 



to ijj 



RG 
K 



It seems that a factor for the discrepancy is 
due to the discretization error despite much improved 
algorithm. Part of the problem could be from the an- 
alytic continuation where very sharp spectral peaks are 
fit with overestimated width. However, it is not clear at 
the moment, given the above values for aqmc^ how such 
discrepancy affects the scaling behaviors. 



u 


Pl/2 


'^1/2/ /3min 




10 


14 


0.135/60 


0.26 


12 


21 


0.091/60 


0.24 


14 


36 


0.048/60 


0.21 



TABLE I: Inverse temperature /3i/2 = l/Ti/j for G{Ti/2, * = 
0) = |Go, bias $1/2 for G(rmin, $1/2) ~ f Go at the minimum 
temperature of Tmin ~ l//3min, and the scaling coefficients 
dnmc derived for U = 10, 12, 14. 



12 



0.4 



0.0 



fa) U = 10, r = 




At= 1/5 
•- fj=m 

0- 0=56 

- 0=56 with Aeq(w) 
^ 13=24 

- fi=24 with Aeq(tj) 
13=] 6 
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(b) No Pade correction 
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FIG. 10: (Color online) Conductance for U — 10 ai inverse 
temperatures /3 — 16, 24, 36 and 60. (a) Conductance curves 
calculated using equilibrium spectral function are shown with 
dashed lines. The narrower zero-bias anomaly width (about 
50 — 60 % from the equilibrium Kondo scale) indicates re- 
duction of the Kondo effect at finite bias, (b) Conductance 
without the Pade correction. The difference is insignificant. 



results are much more robust than spectral function cal- 
culations. 

V. CONCLUSIONS 

Nonequilibrium imaginary-time theory has been for- 
mulated by introducing complex chemical potentials via 
the Matsubara voltage. It has been shovifn that the 
imaginary-time Green's functions upon analytic contin- 
uation are equivalent to the real-time Green's functions. 
For numerical analytic continuations we have given de- 
tailed discussions on the analytic structure of nonequi- 
librium spectral functions and generalized spectral rep- 
resentation in the strongly interacting regime. This for- 
malism has an advantage of having familiar mathemati- 
cal structure as in equilibrium theory and can be readily 
adopted for established equilibrium computational tools 
such as quantum Monte Carlo (QMC) method. 

Application of the Hirsch-Fye QMC method to 
nonequilibrium has produced the nonlinear conductance 
physics where the Kondo resonance is strongly reduced 
by the external bias. The conductance peak has a re- 
duced width from the prediction of equilibrium calcula- 
tions. By correcting discretization errors in QMC, reli- 
able conductance has been obtained as a function of tem- 
perature and bias. Using a scaling form of the conduc- 
tance, we obtained coefficients to the leading temperature 
and bias dependent terms and their ratio aqmc ^ 0.2, 
larger than the perturbative prediction apert — 0.15 
but within the experimental values. This suggests that 
nonperturbative effects lead to more rapid quenching of 
Kondo resonance at finite bias. 

This work shows that the imaginary-time theory pro- 
vides an effective computational tool, along with other 
numerical methods, in the fast-evolving field of nonequi- 
librium quantum many-body theory. This method has 
been applied to complex molecular quantum dot sys- 
tems^ and can be readily extended to bulk nonequilib- 
rium using the dynamical mean-field theory, and quan- 
tum dot systems of complex geometry. 



In a non-interacting resonant model with a rigid spec- 
tral function independent of temperature and bias, ao = 
0.25 can be easily obtained. Small a values can be in- 
terpreted as strong temperature dephasing of the Kondo 
resonance compared to that from bias voltage. Our ra- 
tios aqmc have larger values than the perturbative esti- 
mate^ apert = 0.15 from the effective Fermi liquid ex- 
pansion, aqmc falls within the experimental estimates 
which vary over a range of values, aexp = 0.05^, 0.10"^, 
0.252.. We note that the conductance peak dies away at 
bias much smaller than where the spectral shoulders ap- 
pear in Fig.O Therefore, such fine structures should not 
affect the scaling behavior. In general, the conductance 



VI. ACKNOWLEDGEMENTS 

The author thanks Ryan Heary, Karyn Le Hur, 
Frithjof Anders, Akira Oguri, Reinhold Egger and 
David Goldhaber-Gordon for helpful discussions. Spe- 
cial thanks go to Andreas Dirks and Thomas Pruschke 
who pointed out the discretization problems in QMC 
and made their continuous-time QMC calculations avail- 
able. Author is grateful for financial support from the 
National Science Foundation with the grant numbers 
DMR-0426826, DMR-0907150 and computing resources 
at CCR of SUNY Buflfalo. 



13 




FIG. 11: (a) Fourth order vertex correction to the self-energy with the imaginary-frequency labels, (b) The same diagram in 
the real-time theory. Numerical labels denote scattering states, i.e. 1 = (ai,fci,(7i). (c) Time-ordering for the greater Green's 
function with the interaction times (si, S2) with si extending to f = T on the upper Keldysh branch, (d) The same diagram 
with si on the lower branch. 



Appendix A: Vertex correction and branch cuts for = itpm — $ variable 



To examine the general analytic structure of the spectral representation of the imaginary-time self-energy, we go 
beyond the second-order contribution. Here we discuss that the spectral representation is expressed as 

S(zw„,z(p™) = > / de- — , (Al) 

^ J ILOn - ^{ifm - $) - e 

and how the analytic continuation iipm — > $±177 is taken. Up to the second order, the spectral function cr^(e, iipm, — ^) 
does not have any dependence on iipm — To see how ^^(e, iip,n — *&) acquires the iipm — dependence in the high 
order perturbation, we examine the vertex correction shown in Fig. Illf a). First, we express the polarization diagram 
Po{iiym) as 

Paiivm) = ■^^G'o(ia;„ + w,„)Go(iw„) (A2) 
_ ^ V V (de (de r„.r„J.g(6i)p|g(62)P 



Here we introduce short-hand notations, = / de;, pi = {Ti/'iT)\g{ei)\'^ , h = + ^{ifm ~ Then 
uc \ 1 [ f P^P^ f fPiP^ih-fi) 

P ^ J1J2 (l-^n + ll^ni - ei)(«Wn - £2) Jl J 2 l-^m " Cl + £2 



which can be rewritten as 



7=0, ±1' 

with 

^ J^J^PiP2{.f2 - fi)S{e ~ei + £2)5^,^1-^2. (A6) 

The diagram in Fig. [TIJa) becomes 

E(a)(iw„) = ^ E ^o(*'^" ^ iiym)Go{iujn - iJ^™ - Wp)Po(«i^m)Go(ia;„ - ivp)Po{ivp). (A7) 
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Summing over gives the partial factor 



P1P2A3 



1J2J3 



/1-/2 



/i + 



{ivp + I2- h){i^^n - i^p - £3 - £2) («w„ - £3 - ei)(ia;„ - iVp - £3 - £2) 



(A8) 



with the Bose-Einstein fmiction 



term proportional to /i — /2 with the remaining factors in Eq. (jA7l) . we have 
P „ (^'^p + £2 - £i)(iw„ - ii/p - £3 - £2) 



1) ^ = [e^*^^* — 1] ^. Performing the summation on ivp on the first 



(A9) 



4 J5 £4 ~ £2 — £3 



"5 - ni^2 



£5 + £2 - £1 V tw. 



1 



1 



£1 - £3 



+- 



ns, + /2+3 



ILOn - £1 + £2 - £4 

+ h 



with /i = 1 - /, = (e 



£2 - £3 - £5)(«W„ - £1 - £3) {iUJn - £4 - £5)(i'^n - £l + €2 - £4). 

^ = [e-'3('=i-"i*/2) ^ 1]^^. This expression can be reduced to the form 



(AlO) 



(All) 



by repeatedly using 



lUJn — Zi IUJ„ 



Z2 



Zl - Z2 



IU}„ 



Zl 



Z2 



(A12) 



It can be shown that the above form can be deduced for other types of high order perturbation diagrams. The form 
Eq. (jAllI) could have been anticipated from the Hamiltonian Eq. (jlip where i^m ~ serves as a parameter and the 
equilibrium imaginary-time theory has a similar spectral representation for electron self-energy due to the causality. 

However there remains an important question concerning the direction of analytic continuation of ii^m ^- The 
denominator (za;„ — ^{i(pm — <&) — e)~^ in Eq. (lAlip does not pose a problem regarding the direction i(pm — >■ $ + iO+ 
or iipm — )■ $ — iO^ due to the finite imaginary number in iw„. However, the factor B^{e^i^pm — ^) contains energy 
denominators such as (£4 — £2 — £3)""'^ in Eq. (jAlOp . which may result differently depending on the direction of the 
continuation i^Pm — ^ 4' ± iO'*'. 

To resolve this issue, we examine how such energy denominators behave in the Keldysh real-time theory. We 
consider the same fourth order diagram as shown in Fig. [TlTb). For a specific time ordering of Fig. [TTJc) for t > 0, 
its partial contribution to the self-energy can be expressed as 

° dsi f ds2{d{t)d^ (S2)) {d{s2)d\si)) (dt (O)d(si)) {d{s2)d^ (0)) {d^ {s2)d{0)) {d^ {t)d{si)) {d{t)d^ (si)) (A13) 
T Jo 



rt 

dsi / ds2 
T Ja 



n 



/l/2/3/4/5/6/7e 



-ieiSi+iiS4(si-S2)+ie7(s2-t)-ie3S2+ie2S2+i(iS5-«6)(t-si) 



(A14) 



with the continuum labels defined in Fig. Illf c). Here we take the limit T — >■ — oo as prescribed by Gell-mann and 
Goldberg^^ by taking the T- integral r] dT e^^ . Also, performing the integrals on si and S2 we get 



n [p 



1=1.7 



/1/2/3/4/5/6/7 



g-it(e3-e2+e4-e5-l-e6) _ g-it{ee+e7-e5) 

(£1 - £4 + £5 - £6 + «'7)(£3 - £2 + £4 £7) 



(A15) 



Here, the bias dependence is only in the statistical factor /1/2 ■ ■ ■ fr- A different contribution to is given in 

Fig. ITlT d) with si extending to T on the lower Keldysh branch. Its contribution can be similarly calculated as 



n [p 



i=l,7 



/1/2/3/4/5/6/7 



(£1 - £4 + £5 - £6 + «??)(£3 - £2 + £4 - £7) ' 



(A16) 
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with the negative sign coming from different Wick contraction. The only difference from the previous expression is 
the statistical factor. Regarding the convergence factor irj, we are concerned with the contribution A 



(/1/4/5/6 — fififbfe)- 



i=l,4,5,6 ai 

Within the constraint given by the (5-function, 

/1/4/5/6 - /1/4/5/6 = /1/4/5/6 e'^(^^+^«) 



-it (£3 -£2 +£4 -£5 +£6) 



-it(£6+£7-£5) 



£3 - £2 + £4 - £7 



-(5(ei-e4 + £5-£6)- (A17) 



(04+06)* 



-/3(Qi+a5)<E> 



(A18) 



In equilibrium <J> = 0, A = and the energy integral becomes principal value integral and the presence of ir] becomes 
irrelevant. The same holds in nonequilibrium for single quantum-dot systems by applying the same argument in 
Section H] and Fig. [31 States (15) play the role of incoming state \n) and (46) the outgoing state |m) in Fig. [31[a). For 
example, for (aia^) ~ {LR) and (a^ae) = (LL) there exists a permutation of reservoir labels in the a-summation 
of Eq. (|A17|) . (aias) — {LL) and (aiae) ~ (LR) without permuting the energy variables (ei, €4, €5, eg) and changing 
the factor rir4r5r6. Therefore the expression A becomes zero for nonequilibrium and the integrals of energy poles 
around the real axis due to irj can be replaced by principal integrals. Finally the analytic continuation of iipm — ^ ^ 
can be taken as 



and the subtlety of the analytic continuation (iipm ' 
We write the total imaginary-time self-energy as 



<^> + iO+) + (iipr, 
$) is resolved. 



$-iO+)] 



(A19) 



(A20) 



with the function Q expressed as a Fade quotient 



?-y(e,z) = 



In particle- hole asymmetric limit, one also needs to consider the constant term in addition to Eq. (jA20[) . 

S(lW„, liprn) = S [lipm - ^') + > ^ / «£ — 

7 

^'^{ifiyi — ^) is represented by another Fade approximant as 

l + c(l)z + r-(2)z2 



$) 



(A21) 



(A22) 



E°(z) = 



z' + 



l + rf(l)z + rf(2)z2 + 



(A23) 



In this work, we have only considered the particle-hole symmetric limit and the constant self-energy term TP(iipm — 4") 
has not been included in the analytic continuation. 
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